/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Definitions of the Snapshot class. \ingroup sfrhist

// $Id$

#include <iostream>
#include <list>
#include <set> 
#include <string>
#include "preferences.h"
#include "biniostream.h"
#include "counter.h"
#include "constants.h"
#include "snapshot.h"
#include "message.h"
#include <cmath> 
#include "boost/lambda/lambda.hpp"
#include "boost/lexical_cast.hpp"
#include "mcrx-debug.h"
#include <sys/stat.h>
#include "misc.h"

using namespace std;
using namespace blitz;

vector<string> mcrx::Snapshot::species_names_;

/** Constructor just initializes some variables, it doesn't do any work. */
mcrx::Snapshot::Snapshot(const Preferences& gp, const Preferences& p, 
			 const std::string& file, 
			 const std::string& name) :
  prefs_(gp), user_prefs_(p), name_(name), 
  // indicates units are unknown
  massunit_ (0), lengthunit_ (0), timeunit_ (0), a_(0), time_(0)
{
  // initialize the sets vector (i wish these could be consts but I
  // don't know how to fill the vector using an initializer...)
  sets_.push_back(&gas_particles);
  sets_.push_back(&halo_particles);
  sets_.push_back(&disk_particles);
  sets_.push_back(&bulge_particles);
  sets_.push_back(&nstar_particles);
  sets_.push_back(&bh_particles);

  // initialize the star_species set
  star_species_.insert(disk);
  star_species_.insert(bulge);
  star_species_.insert(nstar);

  // make sure species names are inited
  species_name(gas);

  load(file);
}

// Goes through particles looking for one with the correct ID and
// returns the species and index.
pair<mcrx::Snapshot::species, int> mcrx::Snapshot::find_ID(int id) const
{
  for (int s=0; s<n_species; ++s) {
    try {
      return make_pair(species(s), find_ID(species(s), id));
    }
    catch (particle_set_generic::illegal_id&) {};
  }
  // if we come out here, the ID was not found. Now we throw for real.
  throw particle_set_generic::illegal_id();
}


/** Loads a snap file into the snapshot object.  Data are kept in
    snapfile units, but the unit values are set (if possible) so that
    a subsequent call to makephysicalunits can convert them to
    fiducial units. A simple heuristic is used to determine if the
    file is a GADGET snapshot, or if nbodycod=="pkdgrav", it tries to
    read it as a tipsy file.  */
bool mcrx::Snapshot::load(const string& snapfn)
{
  // see if path is a directory
  struct stat sb;

  if (stat(snapfn.c_str(), &sb) == 0 && S_ISDIR(sb.st_mode)) {
    vector<string> files=word_expand(snapfn+"/*");
    cout << "Snapfile given is a directory. Loading files:\n";
    for(int i=0; i<files.size(); ++i)
      cout << '\t' << files[i] << '\n';
    // load ourselves as the first file, so we get header (first file
    // might be another directory, after all...)
    bool success=load(files[0]);
    // and now concatenate the rest onto us
    for(int i=1; i<files.size(); ++i) {
      Snapshot s(prefs_, user_prefs_, files[i], string("temporary"));
      *this += s;
    }
    return success;
  }

  if (snapfn.find(".hdf5")!=string::npos) {
    // hdf5 file
#if defined (HAVE_LIBHDF5) && (HAVE_LIBHDF5_CPP) && (HAVE_H5CPP_H)
    return loadgadget_hdf5(snapfn);
#else
    cerr << "Not compiled with HDF5 support.";
    return false;
#endif
  }

  if ((snapfn.find(".fits")!=string::npos) && 
      (user_prefs_.getValue("nbodycod",string())=="enzo")) {
    // fits file is converted enzo file (this needs to be made more general)
    return loadenzo(snapfn);
  }
  if(isgadget(snapfn))
    return loadgadget(snapfn);
  else if (user_prefs_.defined("nbodycod") && 
	   (user_prefs_.getValue("nbodycod",string())=="pkdgrav")) {
    return loadtipsy(snapfn);
  }
  else if (user_prefs_.defined("nbodycod") && 
	   (user_prefs_.getValue("nbodycod",string())=="ART")) {
    return loadart(snapfn);
  }
  else {
    cerr << "Don't know what format file " << snapfn << " is..." << endl;
    return false;
  }
}

  
mcrx::Snapshot&
mcrx::Snapshot::operator+=(const Snapshot& rhs)
{
  // add the particle sets
  gas_particles += rhs.gas_particles;
  halo_particles += rhs.halo_particles;
  disk_particles += rhs.disk_particles;
  bulge_particles += rhs.bulge_particles;
  nstar_particles += rhs.nstar_particles;
  bh_particles += rhs.bh_particles;

  // everything else is ASSUMED to be the same, so we are done.
  return *this;
}

  
/** Converts the units of the snapshot to physical units. The unit
    variables are updated, so the function can be called several times
    without anything bad happening. */
bool mcrx::Snapshot::makephysicalunits()
{
  if (massunit_ == 0 ||lengthunit_ == 0 ||timeunit_ == 0)
    // we don't have sufficient information to perform unit conversion
    return false;

  time_ *= timeunit_;

  // The particle sets convert their own units. That's a relief!
  for(int i=0; i<n_species; ++i)
    sets()[i]->convert_units(massunit_, lengthunit_, timeunit_);

  // Update the unit values so they reflect the current units
  // This means that running this routine once more will multiply
  // everything by 1...
      
  lengthunit_ = 1;
  massunit_ = 1;
  timeunit_ = 1;
  units_["length"] = "kpc";
  units_["mass"] = "Msun";
  units_["time"] = "yr";
  units_["temperature"] = "K";

  //value of the gravity const;
  G_ = 4.4967e-24;

  return true;   
}


template <typename T_field, typename T_file>
bool 
mcrx::Snapshot::read_field (binifstream& stream, int n, 
			    vector<T_field>& field,
			    vector<int>& read_order,
			    /// Dummy parameter for template inst
			    T_file,
			    bool check) const
{
  // Resize vector here
  field.resize(n);

  if (n == 0)
    // if the field is empty, there is no padding, so we can't read
    // those or we mess up the file
    return true;

  assert (read_order.size() >= n);
  if (check && !check_field_header(stream, n*sizeof(T_file))) {
    DEBUG(1, cerr << "Error: Field position confusion" << endl;);
    return false;
  }
  for(int i=0;i<n; ++i) {
    T_file temp;
    stream >> temp;

    assert (read_order [i] < field.size());
    field [read_order [i]] = temp;
  }
  if (check && !check_field_header(stream, n*sizeof(T_file))) {
    DEBUG(1, cerr << "Error: Field position confusion" << endl;);
    return false;
  }
  else {
    DEBUG(1, if(check) cerr << "Correctly "; \
	  cerr << "Read field with " << n*sizeof(T_file) << " bytes"<<endl;);
    return true; 
  }
}


// Checks that the header in the stream is consistent with n bytes
bool mcrx::Snapshot::check_field_header(binifstream& stream, int n) const 
{
  int bytes;
  stream >> bytes;
  if (bytes != n) {
    DEBUG(1, cerr << "Error: Unexpected field header. Expected " << n << ", got " << bytes << endl;);
    throw n;
    return false;
  }
  return true;
}


// Skips a field in the snapshot file using the size information in the padding.
bool mcrx::Snapshot::skip_field (binifstream& stream, int expected=-1) const
{
  if (expected==0)
    // if there are no particles, there is no field, so we just return
    return true;

  int bytes, temp;
  stream >> bytes;
  stream.skip(bytes );
  stream >> temp;
  if (bytes != temp)
    cerr << "Header error when skipping snapshot field!" << endl;
  else {
    DEBUG(1, cerr << "Skipped field of " << bytes << " bytes"<< endl;);
  }

  return bytes != temp;
}

/** Dirt cheap heuristic for figuring out if the file is a gadget
    snapshot. */
bool mcrx::Snapshot::isgadget(const string& snapfn)
{
  binifstream snap;
  int next;
  snap.open(snapfn.c_str());
  snap >> next;
  if(!snap){
    // the file name may not be an actual file, so if we can't open
    // the file we conclude it's not a gadget file

    return false;
  }
  snap.close();
  if(next==256)
    byte_swap = false;
  else if (next == 65536)
    // this indicates a snapshot with the opposite byte order... 
    byte_swap = true;
  else 
    return false;
  return true;
}

/** This function simply tries to find the particles of all the parent
    ID numbers to make sure they exist. */
void mcrx::Snapshot::check_parent_ID_integrity() const
{
  try {
    for (int i=0; i< nstar_particles.pid_.size(); ++i)
      find_ID(nstar_particles.pid_[i]);
  }
  catch (...) {
    cerr << "Parent ID integrity check failed for new stars!" << endl;
    throw;
  }

  // Also check parent IDs for gas particles
  try {
    for (int i=0; i< gas_particles.pid_.size(); ++i)
      find_ID(gas_particles.pid_[i]);
  }
  catch (...) {
    cerr << "Parent ID integrity check failed for gas!" << endl;
    throw;
  }
}


void mcrx::Snapshot::setup_ordering(vector<int>& read_order,
				    particle_set_generic& set,
				    map<unsigned int, int> shuffle,
				    int npart)
{
  // then loop over shuffle map (in ID order) to generate the
  // sequential position in ID order and map this to the position in
  // the file.
  read_order.resize(npart);
  int id_pos = 0;
  for(map < unsigned int, int >::iterator i = shuffle.begin();
      i != shuffle.end(); ++i, ++id_pos) {
    const int id = i->first;
    const int file_pos = i->second;
    assert (file_pos < npart);
    assert (id_pos < npart);
    // Fill read order vector
    read_order[file_pos] = id_pos;
    // and fill ID order map
    set.order_[id] = id_pos;
  }
} 

void mcrx::Snapshot::calculate_temperatures(vector<float>& u_ngas,
					    vector<float>& ne,
					    vector<float>& tp)
{
  // We can now calculate temperature from internal energy and
  // electron density

  // calculate mp/k in simulation units, ie K*timeunit^2/lengthunit^2
  // (1.16e20 is mp/k in K*yr^2/kpc^2)
  const double mpk = constants::mp/constants::k*
    pow(lengthunit_*constants::kpc/(timeunit_*constants::yr),2);  

  for(int i = 0; i < gas_particles.n(); ++i) {
    // calculate temp in K from internal energy. TJ says:
    // T=u*(gamma-1)*mu*m_p/k where gamma=5/3, m_p=proton mass,
    // k=Boltzmann's constant, mu*m_p is the mean molecular weight, where
    // mu=(1+4*yhelium)/(1+yhelium+nume), yhelium=(1-0.76)/(4*0.76)
    // and nume is the electron density in the snapshot
    const double yHe=(1 - 0.76)/(4*0.76);
    const double mu=(1 + 4*yHe)/(1 + yHe + ne[i]);
    const double T= u_ngas[i]*mu*(5./3 - 1)*mpk;
    const double Teff= (u_ngas[i] + tp[i])*mu*(5./3-1)*mpk;
    gas_particles.temp_.push_back(T);
    gas_particles.teff_.push_back(Teff);
  }
}
    

/** Loads a GADGET snapshot file. This function does NOT handle comoving snapshots.  */
bool mcrx::Snapshot::loadgadget(const string& snapfn)
{
  using namespace boost::lambda;
 
  cout << "Loading Gadget snapshot type 1 file: " << snapfn << endl;

  // this is now only used locally in this file because it's obsolete
  Gadgetheader snapfile;
  snapfile.loadfrheader(snapfn.c_str());
  
  time_ = snapfile.time; 

  // syntactic convenience. note that we can't use the n(species)
  // function yet, because it relies on the particle data being loaded.
  const vector<int>& npart(snapfile.npart); 
  const int n_tot = accumulate(npart.begin(), npart.end(), 0);

  // check for Springel BH
  if(npart[bh]>0) {
    cerr << "Snapshot file contains " << npart[bh] << " BH particles\n";
  }

  DEBUG(1,snapfile.printheader(););

  // Get units from Gadget parameter file (NOTE: UNITS ARRIVE AS cgs units)
  // length in kpc
  lengthunit_ = (prefs_.getValue("UnitLength_in_cm",double()))/
    (constants::kpc*1e2);
  // mass in Msun
  massunit_ = (prefs_.getValue("UnitMass_in_g",double()))/
    (constants::Msun*1e3);
  // time in yr
  timeunit_ = prefs_.getValue("UnitLength_in_cm",double())/(prefs_.getValue("UnitVelocity_in_cm_per_s",double())*constants::yr);  

  // If "h-inverse_units" keyword is set to true, Gadget snapshots have
  // h^-1 in the units, so we should multiply the length, time, and mass
  // units by h^-1
  if(user_prefs_.defined("h-inverse_units") &&
     user_prefs_.getValue("h-inverse_units", bool())) {
    T_float h0 (snapfile.HubbleParam);
    cout << "  Using h^-1 units with h=" << h0 << endl;
    lengthunit_ /= h0;
    massunit_ /= h0;
    timeunit_ /= h0;
  }

  G_=0; // Note: _G will be incorrect until makephysicalunits() is run, we can figure out a _G here and we should  
  units_["length"] = boost::lexical_cast<string>(lengthunit_)+"*kpc";
  units_["mass"] = boost::lexical_cast<string>(massunit_)+"*Msun";
  units_["time"] = boost::lexical_cast<string>(timeunit_)+"*yr";
  units_["temperature"] = "K";


  // Now that header info is set up, open the file
  binifstream snap;
  snap.open(snapfn.c_str());
  if(!snap) {
    cerr << " Error opening SNAP file!\n";
    return false;
  }
  else
    cerr << "Reading GADGET snapshot file " << snapfn << endl;
  if (byte_swap) {
    cerr << "  swapping byte order" << endl;
    snap.swap();
  }
  
  // Because it's easier to put the particles in the arrays once we
  // have figured out the ordering, we skip ahead in the file to read
  // the ID numbers, then reopen and fill the arrays.
  //snap.skip(264);   // header
  skip_field(snap); // header
  skip_field(snap); // pos
  skip_field(snap); // vel

  // There are 3 different "order" containers:
  // the "order" maps map ID->position in the *vector*, for retrieving by ID
  // the "shuffle" maps map ID->position in the *file*, only used to set up:
  // the "order" *vectors*, which map pos in the file->pos in the vector
  // (the positions are guaranteed to be contigous, so a vector can be used)
  // the first is saved with the particle data, the rest only used in this func

  vector<map<unsigned int, int> > shuffle(n_species);
  vector<vector<int> > read_order(n_species);

  check_field_header(snap,n_tot*sizeof(int));
  for(int s=0; s < n_species; ++s) {

    // Set up shuffle map for species s
    for(int i=0;i<npart[s]; ++i) {
      unsigned int temp;
      snap >> temp;
      shuffle[s][temp] = i; // ID->file pos
    }

    // Check that all ID numbers are unique (at least for a given species)
    if(shuffle[s].size() < npart[s]) {
      cerr << "FATAL: Duplicate ID numbers" << endl;
      throw particle_set_generic::illegal_id();
    }

    setup_ordering(read_order[s], *sets()[s], shuffle[s], npart[s]);
  }
  check_field_header(snap,n_tot*sizeof(int));
    

  // Now that we have the read_order set up, it's easy to read the fields.
  // Reopen file and start over

  snap.close();
  snap.open(snapfn.c_str());
  skip_field(snap); // header

  // Read positions
  check_field_header(snap, 3*n_tot*sizeof(float));
  for(int s=0; s < n_species; ++s)
    read_field (snap, npart[s],sets()[s]->pos_, read_order[s], 
		TinyVector<float, 3>(), false);
  check_field_header(snap, 3*n_tot*sizeof(float));

  // Read velocities
  check_field_header(snap, 3*n_tot*sizeof(float));
  for(int s=0; s < n_species; ++s)
    read_field (snap, npart[s],sets()[s]->vel_, read_order[s], 
		TinyVector<float, 3>(), false);
  check_field_header(snap, 3*n_tot*sizeof(float));

  // Read IDs
  check_field_header(snap, n_tot*sizeof(int));
  for(int s=0; s < n_species; ++s)
    read_field (snap, npart[s],sets()[s]->id_, read_order[s], int(), false);
  check_field_header(snap, n_tot*sizeof(int));

  // Masses are a little more complicated because they might be
  // defined to be constant and not be in the field.

  // We don't bother to check the length of the field because we would
  // have to figure out how many are in the field in advance...
  snap.skip(4 );
  for(int s=0; s < n_species; ++s)
    if(snapfile.mass[s] == 0)
      // Read individual masses
      read_field (snap, npart[s],sets()[s]->m_, read_order[s], float(), false);
    else {
      // Mass is constant, just fill the vector
      sets()[s]->m_.assign(npart[s], snapfile.mass[s]);
      DEBUG(1, cerr << "Constant mass for "<<npart[s]<< " particles"<<endl;);
    }
  snap.skip(4 );


  // Read Internal Energy (gas only)
  // This is used to calculate temp later so for now read it into a temporary
  vector<float> u_ngas;
  read_field (snap, npart[gas], u_ngas, read_order[gas], float());


  // Read densities (currently not saved)
  skip_field (snap, 4*npart[gas]);
  
  // Read Electron Density, if cooling is enabled. This is also used
  // for the temp calculation later so it's also in a temporary.
  vector<float> ne;
  if(snapfile.flag_cooling==1)
    read_field (snap, npart[gas], ne, read_order[gas], float());
    
  
  // Read Neutral Hydrogen Density (gas only, currently not saved)
  if(snapfile.flag_cooling==1) {
    skip_field (snap, 4*npart[gas]);
  }
  
  // Read Hybrid particle data 1 (formed stellar mass).

  // This field should never be present, ie flag_stargens should
  // always be set, because all gadgets are currently using the
  // stochastic SFR.
  if(snapfile.flag_sfr &&!snapfile.flag_stargens) {
    // we should check this at the start
    cerr << "FATAL: This file contains hybrid particles. Can't read." << endl;
    exit(1);
  }
  
  // Read Hybrid particle data 2 (mass of cold clouds, gas only,
  // currently not saved)
  if(snapfile.flag_multiphase==1) {
    skip_field (snap, 4*npart[gas]);
  }

  // Read smoothing lengths
  read_field (snap, npart[gas], gas_particles.h_, read_order[gas], float());
  
  // Read SFR (gas only)
  // Note: This field is in M_Sun/year, unlike the others which are in
  // gadget units
  if(snapfile.flag_sfr==1) {
    read_field (snap, npart[gas], gas_particles.sfr_, read_order[gas], float());
  }
  else
    gas_particles.sfr_.assign(npart[gas], 0.);

  // back-convert SFR to gadget units so that it can get converted in
  // the general setup routine (this is retarded)
  transform (gas_particles.sfr_.begin(),gas_particles.sfr_.end(),
	     gas_particles.sfr_.begin(),
	     _1 * (timeunit_/massunit_));


  // Read Formation Time (new stars only)
  DEBUG(1, cout << "Reading stellar age." << endl;);
  if(!snapfile.flag_stellarage) {
    // this should never happen but there could be no stellar particles at all
    cerr << "Warning: The file contains no stellar ages." << endl;
  }
  else {
    // if we're using the stochastic star formation, this field only
    // applies to the new stars, otherwise it's the mean stellar age
    // of both new stars and gas particles. But we have already
    // checked that we use stochastic SFR.

    // This field is only present if npart[4] > 0
    if(npart[4]>0)
      read_field (snap, npart[nstar], nstar_particles.tf_, 
		 read_order[nstar], float());
  }
     

  // Read feedback reservoir energy (originally called "turbulent pressure"). 

  // This is used to calculate the "effective" temperature by adding
  // to the internal energy
  vector<float> tp;
  if((snapfile.flag_feedback==1) && (!snapfile.is_springel_file)) 
    read_field (snap, npart[gas], tp, read_order[gas], float());
  else
    tp.assign(npart[gas], 0.);

  // If we are reading a Springel file, metallicity is here.
  if(snapfile.is_springel_file && snapfile.flag_metals) {
    DEBUG(1,cerr << "Reading metallicity." << endl;);
    check_field_header(snap, (npart[gas]+npart[nstar])*sizeof(float));
    read_field (snap, npart[gas], gas_particles.z_, 
		read_order[gas], float(), false);
    read_field (snap, npart[nstar], nstar_particles.z_, 
		read_order[nstar], float(), false);
    check_field_header(snap, (npart[gas]+npart[nstar])*sizeof(float));
  }
  
  DEBUG(1,cerr << "Reading potential." << endl;);
  // Gravitational potential (all particles, not used)
  if (snapfile.flag_snapaspot)
    skip_field(snap, 4*ntot());

  // If we are reading a Springel file, BH mass and BH mass accretion
  // rate are here. The units are read in as GADGET units and
  // converted to SUNRISE units.
  if (snapfile.is_springel_file) {

    read_field (snap, npart[bh], bh_particles.mbh_,
		read_order[bh], float());
    DEBUG(1,cerr << "BH masses are (in GADGET units)"	\
	  << setiosflags( ios::scientific );		\
	  for (int i=0; i<npart[bh]; i++)		\
	    cerr << ' ' << bh_particles.mbh_[i];	\
	  cerr << endl;);
  
    read_field (snap, npart[bh], bh_particles.mdotbh_,
		read_order[bh], float());
    DEBUG(1,cerr << "BH mass accretion rates are (in GADGET units)"; \
	  for (int i=0; i<npart[bh]; i++)			     \
	    cerr << ' ' << bh_particles.mdotbh(i);		     \
	  cerr << endl << setiosflags( ios::fixed );); 
  }

  // If we are reading a TJ file, metallicity is here.
  if((!snapfile.is_springel_file) && snapfile.flag_metals) {
    read_field (snap, npart[gas], gas_particles.z_, 
		read_order[gas], float());
    read_field (snap, npart[nstar], nstar_particles.z_, 
		read_order[nstar], float());
  }

  // If we haven't read any metals at all so far, set it to zero so
  // the vector is correctly sized
  if(gas_particles.z_.empty())
    gas_particles.z_.assign(npart[gas], 0.);

  // Skip fields we don't care about - TJ file only
  if (!snapfile.is_springel_file && snapfile.flag_energydetail) {
    skip_field (snap, 4*npart[gas]); // totRadE
    skip_field (snap, 4*npart[gas]); // totShockE
    skip_field (snap, 4*npart[gas]); // totFbE
  }

  // Read ID numbers of the parents (new stars only)
  if(snapfile.flag_parentID) {
    // If this flag is set, it's a TJ file which has parent IDs in a field
    read_field (snap, npart[nstar], nstar_particles.pid_, 
		read_order[nstar], int());

    // Because there was a bug in TJs code to write this, it sometimes
    // has zeros in it, but it only occucs for the particles that have
    // inherited their parent gas particles ids so in that case we can
    // copy the id number itself
    bool crap=false;
    for_each(nstar_particles.pid_.begin(),
	     nstar_particles.pid_.end(),
	     var(crap)|=(_1 == 0));
    if (crap) {
      cerr << "Warning: Zero parent ID bug detected. Using own ID." << endl;
      for (int i=0; i< npart[nstar]; ++i)
	if (nstar_particles.pid_[i] == 0)
	  nstar_particles.pid_[i] = nstar_particles.id_[i];
    }
  }
  else if (snapfile.is_springel_file) {
    // The Springel files have parent IDs encoded in the ID's themselves.
    // This is true of both gas particles that spawn new stars and
    // second generation new stars (those spawned from a gas particle that
    // has already spawned a new star).
    // We get the parent ID by erasing the higher bits
    transform(nstar_particles.id_.begin(),
	      nstar_particles.id_.end(),
	      back_inserter(nstar_particles.pid_),
	      _1 & 0x7FFFFFFF);

    transform(gas_particles.id_.begin(),
              gas_particles.id_.end(),
              back_inserter(gas_particles.pid_),
              _1 & 0x7FFFFFFF);
  }
  else
    cerr << "Warning: File contains no parent ID information\n";

  // Read new star creation mass (if present) and creation size
  if(snapfile.flag_starorig) {
    // This is a TJ file which has this field
    read_field (snap, npart[nstar], nstar_particles.mcr_, 
		read_order[nstar], int());

    // Because there was a bug in TJs code to write this, it sometimes
    // has invalid numbers (mostly zeros, sometimes huge numbers) in
    // it. In that case we just set it to the current mass.
    if( mismatch(nstar_particles.mcr_.begin(),
		 nstar_particles.mcr_.end(),
		 nstar_particles.m_.begin(),
		 ((_1/_2 >= 1.0) && (_1/_2 <= 10.))).first != 
	nstar_particles.mcr_.end()) {
      cerr << "Error: Invalid number detected in creation mass!" << endl;
      cerr << "       Setting creation mass to current mass!" << endl;
      nstar_particles.mcr_.assign(nstar_particles.m_.begin(),
				   nstar_particles.m_.end());
    }

    // The creation size is the smoothing length of the gas particle
    // from which the new star particle was spawned.  NOTE THAT THIS
    // IS OVERWRITTEN AT THE END OF THE FUNCTION, but needs to be read
    // here to get the file fields in order.
    read_field (snap, npart[nstar], nstar_particles.h_, 
		read_order[nstar], int());

    // This also has a similar bug and can be set to zero. Check.
    if( find_if(nstar_particles.h_.begin(),
		nstar_particles.h_.end(),
		((_1/prefs_.getValue("SofteningStars",double()) < 1e-6 ) ||
		 (_1/prefs_.getValue("SofteningStars",double()) > 1e6 ))) !=
	nstar_particles.h_.end()) {
      cerr << "Error: Invalid number detected in OrigHsml!\n"
	   << "       Using zero (but it's unused anyway so don't worry)." 
	   << endl;
    }
  }

  if(!snap && !snap.eof()) {
    cerr << "There were problems reading the file!\n";
    return false;
  }
  if((! snap.eof())&&(snap.skip(1).eof()))  
    ;//cerr << "Correctly Read to EOF" << endl;
  snap.close();

  // This concludes the reading of the file


  calculate_temperatures(u_ngas,ne,tp);

  init_disk_bulge_halo();

  
  return true;
}

void mcrx::Snapshot::assign_softening(species s, const string& name,
				      vector<T_float>& h)
{
  const string softkw("Softening"+name);
  const string maxkw("Softening"+name+"MaxPhys");
  if(n(s)>0) {
    const T_float r= prefs_.defined(softkw) ? prefs_.getValue(softkw,double()):1.0;     
    if(!prefs_.defined(softkw))
      cerr << "Warning: Keyword \"" << softkw
	   << "\" not found, using default value 1.0" << endl;
    else
      cout << "Assigning softening " << r << " to " << name << " particles\n";
	
    h.assign(n(s), r);
	     
    if(prefs_.defined(maxkw) && 
       user_prefs_.getValue("nbodycod",string())=="GADGET") {
      // check that we don't go over the max physical softening if we
      // are using gadget
      const double maxphys=prefs_.getValue(maxkw, double());
      for(int i=0; i<h.size(); ++i)
	h[i] = std::min(h[i], maxphys/a_);
      }
  }
}
  
void mcrx::Snapshot::init_disk_bulge_halo()
{
  // There are no specifications for softening lengths for the
  // disk/bulge/halo particles in the snapshot. They are in the
  // parameter file, so if those keywords are defined in the
  // preferences, we use those values.
  // Because of the persistent problems with using the creation
  // smoothing length for the nstar particles, let's just initialize
  // it to the softening length right here.
  assign_softening(disk, "Disk", disk_particles.h_);
  assign_softening(bulge, "Bulge", bulge_particles.h_);
  assign_softening(nstar, "Stars", nstar_particles.h_);
  assign_softening(bh, "Bndry", bh_particles.h_);
  assign_softening(halo, "Halo", halo_particles.h_);

  // init the metallicities and formation times of the disk/bulge pops
  // to zero and Nan, respectively (if they exist in the file, they are overwritten)
  disk_particles.z_.assign(n(disk), 0.);
  disk_particles.tf_.assign(n(disk), blitz::quiet_NaN(T_float()));
  bulge_particles.z_.assign(n(bulge), 0.);
  bulge_particles.tf_.assign(n(bulge), blitz::quiet_NaN(T_float()));

  // init the parent id numbers of the disk/bulge pops to own ids.
  disk_particles.pid_.assign(disk_particles.id_.begin(),
			     disk_particles.id_.end());
  bulge_particles.pid_.assign(bulge_particles.id_.begin(),
			      bulge_particles.id_.end());

  // init the disk/bulge creation mass to current mass 
  // REMEMBER to deal with the mass loss issue somewhere
  disk_particles.mcr_.assign(disk_particles.m_.begin(),
			     disk_particles.m_.end());
  bulge_particles.mcr_.assign(bulge_particles.m_.begin(),
			      bulge_particles.m_.end());
  if(nstar_particles.mcr_.empty()) {
    cout << "No creation mass for star particles in file -- assuming current mass" << endl;
    nstar_particles.mcr_.assign(nstar_particles.m_.begin(),
				nstar_particles.m_.end());
  }
}


string mcrx::Snapshot::species_name(species s)
{
  // Check if name vector is initialized, if not do.
  if(species_names_.empty()) {
    species_names_.push_back("gas");
    species_names_.push_back("dm");
    species_names_.push_back("disk");
    species_names_.push_back("bulge");
    species_names_.push_back("nstar");
    species_names_.push_back("bh");
  }
  return species_names_[s];
}

